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Abstract 

We use the known values of the two loop Wilson coefficients and the three 
loop anomalous dimension matrix 7(n) to perform a next-to-next-to leading order 
(NNLO) calculation of ep deep inelastic scattering. Because 7(n) is only known 
for a few values of n, the method of average reconstruction has to be used, which 
leaves 102 effective experimental points, for 12 parameters: the QCD mass A, 
and 11 initial values for the moments of the structure functions. The data points 
spread in the range of momenta 2.5 GeV^ < Q'^ < 230 GeV^. 

The x^/d.o.f. decreases substantially when going from LO to NLO, and 
also from NLO to NNLO (although only a little now) to /d.o.f. = 79.2/(102 - 
12). The favoured value of A is 

A{nf = 4, 3 loop) = 282.7 ± 35.1 MeV, 

corresponding to the value of the coupling at the Z mass of 

^(3 ioop)|,^2^ ^ Q^^^2 ± 0.0024. 

The calculation, which constitutes a very precise test of QCD, includes target 
mass corrections; the error takes into account experimental errors and higher twist 
effects among other estimated theoretical errors. 
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1 Introduction 

Deep inelastic scattering (DIS), in particular of electrons (or muons) on protons, constituted one of the first 
probes of hadron structure. The calculation of QCD-induced scaling violations in the structure functionsl^l 
yielded some of the earliest qualitative checks of the quark-gluon theory of hadron interactions, as well 
as providing the first two loop determinations of the strong coupling constant '^1. In this last ref. 2, the 
two loop results for the nonsinglet anomalous dimensions'^'*^' were used, together with the one loop Wilson 
coefficients [''1 to get a calculation of the nonsinglet structure functions to next to leading order (NLO).-^ 
Because both anomalous dimensions and coefficient functions were known completely, it was possible to 
use a method devised by Gross'^l to LO, generalized in ref. 2 to NLO, to obtain an exact, pointlike (in x) 
reconstruction of the nonsinglet part of F2ix, Q^), and of 0:^3(2;, Q^). 

Since then the calculations have been extended to include singlet anomalous dimensions (ref. 6; sec 
ref. 7 for a collection of two loop formulas). Furthermore, the two loop Wilson coefficients are also known^*! 
so that, to perform a NNLO calculation, only the three loop anomalous dimensions are needed. There 
does not exist at present a full calculation of these; but a partial computation has been made available 
recently. [®1 This has been used in some papers (see e.g. ref. 10 and work quoted there) to perform an 
evaluation of deep inelastic structure functions to next-to-ncxt-to- leading order, NNLO. These evaluations, 
however, present a number of shortcomings. First of all, the values of the moments which are not known are 
found by interpolation of the known moments; and likewise the experimental values of the moments, with 
which to compare the theory, are obtained also interpolating and extrapolating experimental data. These 
procedures are dangerous in that they imply unknown systematic errors, thus putting in jeopardy an already 
delicate calculation (NNLO effects are in themselves small, generally speaking). Perhaps more important, 
only nonsinglet structure functions (specifically, xF-^ for ;/ scattering in ref. 10) were considered: while the 
best data exist for the ep (or np) structure function F2, with a strong singlet component. 

In the present paper we improve on this in the following ways. First of all, we extend the calculation 
to include singlet as well as nonsinglet structure functions. This allows us to use both the very precise, 
comparatively low energy, electroproduction SLAG data, and the recent, very high energy HERA results, 
plus some intermediate energy muon-production data. This provides a huge range of values. Secondly, we 
employ the method of Bernstein polynomials. This method was developed in ref. 11 for calculating averages 
of the structure functions around values of x where experimental data are available, in terms of only a finite 
set of known moments: thus at the same time avoiding the problem of interpolating the theoretical values of 
the anomalous dimension, and greatly diminishing the errors inherent to the calculation of the experimental 
input, as only known sets of data points are essentially relevant. 

In this way we find a precise comparison of electroproduction DIS data with theory, obtaining in 
particular a very accurate evaluation of the QCD coupling and mass. With respect to the last we find 

^(3 ioop)(-j\^2-) ^Q ji72 ± 0.0017 (statistical) ± 0.0017 (systematic) 

=0.1172 ±0.0024, 

to which corresponds the value and error of the mass parameter 

A{nf = 4, 3 loop) = 283 ± 35 MeV . 

^ To leading order (LO) one requires the one loop anomalous dimension and the tree level Wilson coefficients. 
Likewise, a NLO calculation uses one loop Wilson coefficients and two loop anomalous dimension and the NNLO 
one (which is the subject of the present paper) implies two loop Wilson coefficients and three loop anomalous 
dimensions. 
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2 The evolution equations 

To describe the calculation, let us establish some notation first. We split F2 into a singlet and a nonsinglet 
part, 

F2{x, Q^) = Fsix, Q^) + Fnsix, Q^). (2.1) 

We also define the gluon structure function, Fg{x, Q^) = xG{x, Q^), G being the gluon density. Finally, we 
form the vector with components Fs, Fq: 



We obtain the moments by projecting these with powers of x: 



I^Ns{n;Q^)= [ dxx''-^FNs{x,Q^); 
Jo 

Jo 



(2.2) 



dxa;"-2Fi(x,Q2); 



S, G. 



If wc write a light cone expansion for the moments, this is given in terms of the Wilson coefficients and 
the anomalous dimensions of appropriate sets of operators (for details on this, see any text on QCD, e.g., 
ref. 12). Let t = \ \og{Q'^ / v"^) , with v the renormalization scale, and a = as{Q^)/4TT, a' = as (Q")/47r. For 
the nonsinglet component, we have 



^iNs{n■, Q"^) = Cns {n,as{Q'^)/'i7T) c^p - dt' jNsia') 

. Jo 

whereas the singlet equations are of a matrix character, 

ft(n; Q2) = J2 a.(Q')/4'r) jrexp - ^ dt' 7(0') 



(2.3a) 



(2.3b) 



The jNS, 7 are the NS, S anomalous dimension and anomalous dimension matrix, and the A^s, A unknown 
matrix elements of certain operators. What the calculations of Larin et al.I^l provide is the expression of the 
gammas to three loops (third order in as) for n = 2, 4, 6, 8.^ 

Before discussing how to use this information, let us write the equations that follow from Eqs. (2.3). 
We will consider in some detail only the singlet component; the NS one is fairly trivial. We start by 
completing a two by two matrix C from Cs, Cq, 



C = 



Cs C12 
Cg C22 



Cii C12 
C21 C22 



by requiring it to commute with the anomalous dimension matrix: 

[C,7]=0 

(the index n will be omitted in some of the formulas to lighten the notation). We define 

7(a) 



D{a) 



ma) 



(2.4) 



^ Note that the conventions of ref. 9 differ from ours here in that their 7's are a factor of two larger than ours: 

c\ ref . 9 / \ here / \ 

27 (n) =7 (n). 
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and write the series expansions, 



C{a) =1 + C«a + C(2) a^ + ... 

7(a) =7(0) a + 7^'^ a^+j^^^ + . . . 
-/3(a) =(3oa^ + pia^ + /?2a* + . . . , 
D{a) =-£>(«) +£>W+£)(2)o+... 



(2.5) 



-1 



£)(0) = ^^(0); = Z± { ^(1) _ ^^(0) 



2A 



2/?o 



7^ 



(2) _ ^^^(1) I f ^± _ ^ ^ ^(0) 



/3o 



2/3o 

/3? /?2 



/3i 



-7' 



7' 



(2.6) 



The commutativity of C and 7, order by order in perturbation theory, imphes the equations 



a 



(0) 

(1) _721 ^(1) 
21 -— "-^12 > 

7l2 



a 



(1) 
22 



(0) (0) 
.(1) , 722 -7ii ^(1). 

7l2 



-c ' + 



•'12 1 



a 



(2.7) 



(0) (0) (1) (1) ^(1) ^(1) 

^(2) _^(2) 722 -7ll ^(2) 722 " 7ll ^(1) , C^ll " ^22 .,(1) 



7l? 



7i? 



7i? 



These equations allow us to find the C21, C22 in terms of the Cn, C12 and the -jij which in turn are taken 
from ref. 9. 

For the moments equations, we define the matrix S that diagonalizes the LO anomalous dimension 
matrix, so that 



We take it to be 



D 



(0) 



d- — d+ 



, Sii = det S = 1, 



d--d+ j 



and we define also 



etc. 



Here the matrix M is defined to be such that 

^{a^'°' M(a)C-i(a)Ma)}=0, 

and it has the series expansion 

M{a) = 1 + Af(i)a + Af^^^o^ + • • • . 
From Eqs. (2.8c), (2.9) we find the explicit values of the M: 



(2.8a) 



(2.8b) 



(2.8c) 



(2.9) 



D 



(1) 
11 



V 1 + (i_ - (i+ 



(1) 
21 



(1) 
22 



(2.10a) 
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d+ — d- 



5g) + (M(i)5(i))2i 



2 + rf_ - rf+ 

The evolution equation may then be written as, to NNLO, 



+ (M(l)5(l))22 



/x(n, a) =C{n, a)S{n)M{n, a)"^ {ao/a)'^ '"^ M(n, ao)6'"^(n)C(n, ao)" V(?^, ao) 
a =as (Q^ ) /47r, ao = (Qo ) /^ti" 
and as is to be calculated to three loops: 



(2.10b) 



(2.11) 



with 



An 



L=lo!. 



0,\ogL , pflog^L -pflogL + p2f3o- 



91 



11 



Pi -1L)2 - — n/, 132 - ^ + — 



325„2 



I^Ns{n,Q'^) 



-7<,°J;(n)/2/3o 



We work in the MS scheme throughout this paper. For the NS component the corresponding equation is 

V".s(Q2) 

1 + BL'^a,{Q^)/47r + B^r?\as{Q^)/An) 



(2.12a) 



1 + si'^a«(Q2)/47r + (a«(g^)/47r) 



(2)/ 



-UNs{n,Qo)- 



The explicit values of the Bn \ expressible in terms of the 7^5, C^s('^), niay be found in ref. 10. We can 
rewrite (2.12a) expanding the denominator as 



(2.12b) 



x{l + B(^^as{Q^)/A7r - B^^^asiQ^/An 

^ (i?(i)a,(g2)/47r) (B(i)a,(Q2)/4^)) ^^^sin,Ql). 

The difference between (2.12b) and (2.12a) is of order a^, and may be used to estimate the effect of higher 
order (NNNLO) corrections. 

A last preliminary point to be discussed is that of the number of flavours. We will be working 
with the moments /i(ri, Q^), so we have only one momentum variable. We then split the range into the 
following two intervals: 

Q2 < ml (I); ml < (II). 

The values of wc will be using will be much less than , and larger than or of the order of so wc 
need not consider c, t quark thresholds. Then, in region (I) we take n/ = 4 and in region (II), nf = 5. The 
matching will be carried over following the prescription of ref. 13: to NNLO, 

Po'^' log ^^^(^ = iPo'^' - P7)Lh + ^AT^O + S^NLO 

where 



(2.13a) 



5nlo ={b?^ - hr)\ogLn - 6^^ log 



07^ 
P7 ' 



Snnlo =- 



Po'Lh 



(6^+' - 6^)6^ logL, + K+'r - + b7 - b7^' + h 



(2.13b) 



-4- 



-CALCULATION OF ELECTROPRODUCTION TO NNLO AND PRECISION DETERMINATION OF a^- 

Here, 

Lh = log [m2(n/ + 1)/A\nf)] , hi = A//3o 
and m(n/ + 1) is the pole mass of the (n/ + l)th quark. 

3. Experimental input 



3.1 The method of Bernstein polynomials 

Because, for a given value of ^ only a limited number of experimental points, covering a partial range of 
values of x, are available, one cannot simply use the moments equations. A method devised to deal with a 
situation like the present one is that of averages with the (modified) Bernstein polynomials.^ We define the 
(Bernstein) polynomials as 



Pnk{x) 



2r(n+ I) 



„2fe 



r(A; + i)r(n-fc + l)' 



[l-x^] 



,2\n-k 



2{n-k)\r{n- 



n—k 



'r{k + \)r{n - + 1) § i\{n -k- 1)\ 



(-1)' 



(3.1) 



These polynomials have a number of useful properties. First, they are positive and have a single maximum 
located at 

r(fc + i)r(n + |) 
r(fc + i)r(n + 2)' 

they are concentrated around this point, with a spread of 



^ n - 



r(fc + i)r(n + |) 



n 2 



r{k + \)r{n + 2) 



and they are normalized to unity, dxpnk{x) = 1. Therefore, the integral 

Axpnk{x)(p{x) 

represents an average of the function ip{x) in the region Xnk — \/S.Xnk ^ x < Xnk + ^^Xnk', the values of the 
function if{x) outside this interval contribute little to the integral, as Pnk{x) decreases to zero very quickly 
there. Finally, and using the binomial expansion in Eq. (3.1), it follows that the averages with the Pnk of a 
function can be obtained in terms of its even moments: 



Jo 



1 2(n-fc)!r(n 

dxpnk{x)(p{x) = 



i — k 



{-ly 



r{k + i)r(n - fc + 1) § l\{n - k - l)\ '^'*^+''' 

V2i = / dxx'^'-ip{x). 
Jo 

We will thus consider our experimental input to be given by averages 

Fif^\Q') = j\xpnk{x)Ft''\x,Q% 
Jo 

with F2^^\x,Q'^) the experimental structure function. 



(3.2) 



^ Wc call the polynomials modified because, since only even moments are known, we have to consider polynomials 
in the variable . For more details on the method, see refs. 11, 14. 
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3.2 Calculation of the experimental averages 

Even with the method of Bernstein polynomials, the calculation of an average such as (3.2) requires some 
interpolation and extrapolation of F2'^^^\x, Q^). To do so we have used two different methods. In both we 
separate F2 into a singlet and a nonsinglet part, writing 

F2 = Fs + Fns- 

In the first method we use, for each value of independently, a phenomenological expression for the F, 











x)". 










=8,75 GeV. 


99/02/23 


12.40 


M 








F,(NS) 
F.(S) 


















1 













-4 


10-^ 


10-^ 


' 





Figure 1 Phenomenological fit for = 8.75 GeV (logarithmic scale). 
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Figure 2 Phenomenological fit for = 150 GeV (logarithmic scale). 
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A theoretical justification of (3.3) may be found in ref. 7 and work quoted there; but we want to 

emphasize that, in the present paper, (3.3) is to be considered as only a convenient interpolation of the data 
to allow calculation of the integrals (3.2). In fact, the parameters A, B, C, v, ji are taken to be totally free, 
and assumed to be uncorrelated for different values of Q^. Thus, no theoretical bias is induced in the 
dependence. The fits, for a couple of representative values, are given in Figs. 1,2. In Fig. 3 we also show 
the polynomial P3i{x), superimposed on the fit and the experimental values of F2{x, = 8.75 GeV^), as a 
representative case. 

The second method we consider is to use the parameterization of data (that we denote by MRST98) 
given in ref. 15, which includes some theoretical input. Then the integral involved in the Bernstein average 
is evaluated with the help of the parametric expression, for each value of Q^. 

We consider the first method to be the cleanest one; the results found using the MRST98 param- 
eterization are presented mostly to show the insensitivity of the evaluation to the method of obtaining the 
"experimental" averages f'^-^'^\Q'^) ■ 
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Figure 3 Phenomenological fit for = 8.75 GeV, and Bernstein polynomial pai . 

Only tfiis polynomial is shown as an example, although for this value of Qq other Bernstein 
averages are also meaningful. The points at very small x are not shown (compare with 
Fig. 2). 



3.3 Experimental input 

To begin with, we have the limitation that, since the 7s in the QCD evolution equations are only known for 
n = 2, 4, 6, 8, we will be only able to use, for the theoretical calculations, the ten Bernstein polynomials* 

Pnk{x); fc =0, n = 0, 1, 2, 3; 

k =1, n = 1, 2, 3; , , 

3.4 

A;=2, n = 2, 3; ^ ' 

k=3, n = 3. 



"* Actually, (n) is also known for n = 10 but we will not be able to use this because in F2 for electroproduction 
singlet and nonsinglet are not separated. 
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However, we will not have experimental information, for a given Q^, on all of the corresponding averages. 
Thus, for example, at very large values of Q^, the HERA measurements have been carried over only for 
small values of x, so here only the averages with pso and at most p2o can be used. Also, there are generally 
speaking no data for x near unity, so averages with large n — k cannot be used either. As a matter of fact, 
we have only used the six polynomials 

Pio, P20, Pso; 

P2I, P31, PZ2- 

When we take this into account, it follows that we can calculate reliably only a total of 102 Bernstein 
averages, using the experimental data from SLACti^l, BCDMSt^"^!, E665[i*l and HERAli^l. The experimental 
averages, as defined in Eq. (3.2) with the interpolation given by (3.3), are shown in Fig. 4 (together with 
the NNLO fit, see below); the errors there are the errors in the integrals (3.2) induced by the errors in 
the parameters A, B, C, /i; and these errors are obtained, in turn, from the experimental errors of the 
fitted F2. These errors of f!^^'^'' are what we consider the statistical, experimental errors of our input. 
There is however one more point that has to be discussed in this respect. We have, for each Q^, theoretical 
information on four moments, but we use six Bernstein averages: so we are double counting part of our 
information. This is repaired as follows: we consider the values of Q"^ for which we can use more than 
four Bernstein averages, and consider the number above four of these to be "duplicate" information. The 
total number of these duplicates is of 34. Then we renormalize the experimental error by multiplying it by 
^102/(102-34) = 1.22. 



4. Theoretical calculation. Error analysis 

The QCD equations (2.11, 12) only give the evolution of the moments. Therefore, we have to take their 
initial values, at a fixed Q'q 

lJ-s{n,Ql), iiG{n,Ql); liNs{n,Ql)\ n = 2, 4, 6, 8, 

as parameters. For the actual calculations we chose Qq = 8.75 GeV^ because, for this value, there are 
experimental data in a wide range of values of x so the input values are directly constrained by experiment. 
We have considered varying the value of Qq; this provides an indication of the errors induced by the approx- 
imations made for, in an exact calculation, the result should be independent of Qq. Indeed, the variation is 
small; for example, if we take Qq = 12 GeV^ the QCD mass A varies by 11 MeV. 

Of these twelve moments only eleven are free. In fact, because of the momentum sum rule it follows 

that 

1 

lim [m5(2, Q^) + Mg(2, Q^)] = {QD = - E ^9' 

■' g— 1 

where (Qq) is the average quark charge; for n/ = 4, (Q^) = ^; for ri/ = 5, (Q^) = ^. To implement this, 
we may use the equations (2.11) with as(Q^) 0, and that rf+(2) = 0, d_(2) < together with the values 
of the Sij (2) to get the constraint 

16+^ ^ (M((2,ao)5-^(2)C-^(2,ao))yM,(2,Q2) (4.1) 



and the quantities here are to be calculated for rif = 4. So we may use (4.1) to eliminate /iG(2, Qq) in favour 
of /is(2, Qq), leaving 11 free moments, as stated. To these we have to add a further parameter viz., the QCD 
mass A, to a total of 12 free parameters for the fit. 

The procedure is now rather straightforward. We calculate the //(n, Q^) in terms of the /i(n, Qq) using 
(2.11, 12); and once these /i(n, Q^) known, we evaluate the theoretical Bernstein moments FkniQ^) 
Q^, in terms of our twelve parameters. We are then ready to fit the experimental quantities obtained before, 
Pnk'^\Q^)- a fit which is performed with the MINUIT program[^°l. Before presenting the results, however, 
we have to discuss a number of theoretical questions, including refinements and estimates of theoretical 
errors. 
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4.1. Target mass corrections 

Target mass corrections (TMC) are of order nnip/ Q^, for the nth moment (mp being the mass of the proton) 
and must be included in a precision calculation. The corrections can be evaluated exactlyl^^^, to LO in Ug, 
but we here will only take them into account to first order in rn^/Q'^. The reasons are that higher twist 
(HT) corrections, of order nA^ jQ"^^ are not known; and neither are the NLO (in a^) TMCs. Although it 
may be argued that HT are suppressed with respect to TMC by l/A'^c, with Nc = 3 the number of colours in 
the large Nc limit, it is useless to include effects 0{mp/Q^) while the ones of order A'^/Q'^ and 0{asrn^/Q'^) 
corrections are not known. Anyway, we will give results both with and without TMC; the variation is slight, 
and the error, as just estimated, negligible. 
We write (see refs. 12, 21) 

l^lTin, Q') = ms{n, Q') + !^t^!^^^^sin + 2, Q^). (4.2) 

It should be noted that the TMCs for the singlet are negligible compared to those of the NS, which is why 
we only take into account the last. Another point is that (4.2) involves one further moment, /ijvs(10, Q^). 
Since we do not want to add a new parameter to the fit, what we do is to estimate fj,Ns{W,Q^) from the 
fitted expression for F2, Eq. (3.3). Since this is only a correction, it does not matter much if it is not very 

accurate. 

4.2. Correlation between gluon structure function and A 

As is known, in other fits there is a strong correlation between the value of the QCD parameter A and the 
gluon structure function, which leads to instabilities in the values of the last; this is due to the fact that the 

Ha are not directly related to a measured quantity. The problem is however not serious in our case. This is 
because, in our evaluations, xG{x, Q^) is only represented by its three input moments 

Mg(4, Qo), Mg(6, Qo), mg(8, Qo); 

recall that /xg(2, Qo) can be eliminated in favour of /Us(2, Qo) using the momentum sum rule. We have 
also found that the instabilities disappear completely if we simply make the requirement that /LtG(?^,<9^) > 
Mg('^ + 2, Q^), a requirement that follows from the positivity of G{x, Q^) as a probability density for ghions 
with momentum fraction x. In fact, from the positivity of G we can impose other inequalities;!^^] for example, 
considering only linear ones, we must have 

/■I 2(n-kVT(n+-) ""'^ f-lV 

I <^^Pr..i.)Foix, Ql) = r(,Vi)r(n_4i) E nin-k-D^ + + ^' ^ 

We have checked that his is verified by the central values for the /kg in our fits. 

If we allowed the /ig to vary freely this would result in the appearance of a spurious minimum for 
which the lowest moment would be smaller than the higher ones. Nevertheless, even in this case the value 
of A does not vary too much: it only decreases by some 10 MeV. In our fits we of course require positivity. 

4.3. Quark mass dependence 

Since our range goes through the b quark mass threshold, we have to worry about that region. In our 
main calculation we have simply used the matching conditions (2.13), fixing the mass of the b quark from 
the recent determination in ref. 22, correct to 0(af ), 

mh = 5 001^66^ MeV. 

It is to be noted, however, that the use of (2.13) is not enough to take into account exactly the mass 
dependence. To do so, one would have to rewrite (and solve) the renormalization group equations taking 
into account the finiteness of rrif,, a very hard task quite outside our scope here. The fact that the use 
of (2.13) gives only an approximation to the mass dependence results in a discontinuity of the theoretical 
calculation, at NNLO, near = ml, which can be seen clearly in Fig. 4. Since the discontinuity is of 
^ 1%, and afi'ects only the fit to very few experimental points, we consider the use of (2.13) to be the best 
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available procedure. Alternatively, we could simply avoid the region around the b quark threshold and thus 
cut off the interval^ 

ml<Q^ < Ami 

Actually, we remove from the fit the values of the moments for between 20 and 90 GeV^. Again, we 
take the difference between the results with both methods (matching or cutting the threshold) as part of the 
measure of the error due to the b quark mass effect. The rest of the error is obtained varying mb inside its 

error bars, as above. 

4.4. Higher order corrections 

Higher order corrections are of two types. First, we have higher order in A^/Q^, viz., higher twist effects 
(HT). Since so little is known about these we simply take into account phenomenologically the more important 
ones by adding, to UNsi^jQ^), the correction 

t^Nsin, Q^) = n{aAyQ^)nNs{n, Q^) (4.3) 

with a an unknown parameter, to bo fitted, and expected to be of order miity; in our fit we find a ^ —0.20, 
see below. The difference between the results with this new contribution, and the one without, will be the 
estimated error induced by HT effects. Note that our central values (Table 1 below) are obtained without 
including the HT term. This is because, although the value we get for a is reasonable, the fact remains that 
the expression (4.3) used for HT is little more than educated guesswork, and we want to avoid as much as 
possible to introduce biases in our calculation of yl, Ug. 

The second type of higher order corrections will be corrections of relative order a^, i.e., JVNNLO 
corrections. We estimate these by the different results obtained performing our calculations either with the 
fractional expression (2.12b), or with the expanded denominator one, (2.12c). This simulates, at least in the 
expected behaviour at large n, the largest NNNLO corrections; in fact, the largest corrections of order j to 
F2 may be argued to behave as ai 105:^ (1 — a;) and stem from the C^s(n), 7^5 (n) at large n; see for example 
ref. 12, Sect. 4.9ii and references therein. By using our procedure, we are estimating e.g., 

4.5. Dependence on parameterizations 

The dependence of our input data for the Bernstein averages on the parameterizations used to evaluate the 
corresponding integrals is estimated as follows: we calculate using our interpolation formulas, Eq. (3.3), and 
then repeat the calculation with the MRST98 parameterization. The difference (10 MeV for A) will be our 
estimated error due to this source. 

^ We could also consider increasing the minimum value of to be well clear of the c quark threshold, but we have 
not bothered to do so as the experimental points in the corresponding region are not significant; of. Fig. 4. 
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99/04/08 19.56 

NNLO WITH TMC 

X2 = 79.2/(102-12) = 0.88 
aJM = OJ^ 63±0.0016 (exp) 




I 10 10^ 



Figure 4 Experimental values for the Bernstein averages, and NNLO theoretical 
curves, with the gluon structure function constrained and TMC taken into account. 



5. Results. Discussion 

We present in Table 1 a compilation of the results obtained with our calculations at LO, NLO and NNLO, 
with TMCs taken into account; the fit to the data itself is shown, for the NNLO calculation (with TMC) 
in Fig. 4. In this Tabic 1, only statistical errors arc shown; systematic (theoretical) errors will be discussed 
below. We note that the errors given in Table 1 are "renormalized" to take into account the effective number 
of independent experimental points, as discussed at the end of Sect. 3. That is to say, the error presented 
is obtained from the raw error by multiplying it by 1.22. In calculating the x^/d.o.f., however, we have 
considered all the input experimental points, as both the total chi-squared and the number of degrees of 
freedom change by the same amount. 



Order 


A{nf = 4) 


aa(M|) 




LO 


215 ± 73 


0.130 ± 0.005 


212/(102 - 12) 


NLO 


282 ± 40 


0.116 ± 0.0034 


80.0/(102- 12) 


NNLO 


283 ± 25 


0.1172 ±0.0017 


79.2/(102- 12) 



Table 1 



The NLO corrections are very clearly seen in the fit: the x^/d.o.f. decreases from a largish value 
of ~ 2.4 to a very good ~ 0.89 when including these. The fit is so good already at this order that there 
is very little room for improvement when going to NNLO; nevertheless, an improvement is seen. Not in 
the x^/d.o.f., which only decreases minutely (to /d.o.f . ^ 0.88) but certainly in that including NNLO 
corrections leads to a noticeable gain both in the quality of the determination of the coupling, and in the 
stability of the fits. 

The estimated systematic errors, originating from various sources, are shown for the NNLO case in 

Table 2: 



- 11 - 



-J. SANTIAGO AND F. J. YNDURAIN- 



Source of error 


A(nf = 4; 3 loop) 


AA{nf = 4;3 loop) 




No TMG 


292 


9 


0.0006 


Interpol. (MRST98) 


273 


10 


0.0007 


HT 


292 


9 


0.0006 


Quark mass effect 


299 


16 


0.0011 


Q§ to 12 GeV^ 


294 


11 


0.0008 


AfNNLO 


289 


6 


0.0004 



Table 2 



Let us comment on the meaning of the different entries. No TMC means that we have not taken target mass 
corrections into account. The corresponding error is not included when evaluating the overall theoretical 
error because, since we take into account TMC in our central value, the error would be of order TMC^, or 
as X TMC, quite negligible. MRST98 means that we have used the MRST98 interpolation (Subsect. 3.2) 
to calculate the integrals (3.2). HT means that we have taken into account the presence of higher twist by 
adding a contribution like (4.2). The fitted value of the phenomenological parameter a is a = — 0.202±0.030.^ 
"Quark mass effect" means that we have cut off the h quark threshold as discussed in Subsect. 4.3; the error 
in Table 2 takes into account also the variations due to the error in the rah mass. Qq to 12 GeV^ means that 
we take the input moments defined at this value of the momentum, [ii(n,Q\ = 12 GeV^), i = S, G, NS. 
Finally, JVNNLO means that we have fitted with the expanded formula (2.12c) instead of the fractional one 
(2.12b). Note that we only take this into account for the nonsinglet. This is because for large n, which is 
when the higher order corrections are larger, the NS piece of the structure function is the dominating one. 

Composing quadratically systematic (theoretical) and statistical (experimental) errors we find the 
best result for the QCD coupling, 

A{nf = 4, 3 loop) =282.7 ± 25.1 (stat.) ± 24.5 (syst.) = 282.7 ± 35.1 MeV; 

a(3'°°P)(Mz) =0.1172 ± 0.0017 (stat.) ± 0.0017 (syst.) = 0.1172 ± 0.0024; ^^'^"^ 

the corresponding central value for A{nf = 5, 3 loop) is of 0.200 MeV. It is to be noted that composing the 
"theoretical" errors as if they were independent leads to a certain amount of double-counting. Thus, the 
results should be independent of the value of Qq if the calculation was to all orders, so the two last errors in 
Table 2 are connected. We have, however, preferred to play it safe, particularly because these errors are in 
some cases little more than rough estimates. 

Wc next show a table comparing our results to previous determinations for Q!s(M§); these arc taken 
from the review by Bethke^^"^!, including only the processes where the theoretical calculation has been pushed 
to the NNLO level and the experimental data are good enough to make the analysis meaningful.^ We also 
incorporate in the table the recent results in ref. 10: 



® This is a very reasonable value that makes the mass associated with higher twists, Mht = |o| ~ 0.13 GeV of 
the order of typical HT parameters: intrinsic average transverse momentum of the quarks, or inverse proton radius; 
see e.g. the paper of De Riijula et al., ref. 1. 

This excludes e+e" annihilations. 
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Process 


Average 

or (~)^ rPinp'P frrpVl^ 


a«(M|) 


DIS; V, Bj 


2.5 


0.122 


DIS; V, GLS 


3 


0.115 


r decays 


(1.777)2 


0.119 


Z —>■ hadrons 


(91.2)2 


0.124 


DIS; u, xFa 


5 - 100 


0.117 ± 0.010 


our result 


2.5 - 230 


0.1172 ±0.0024 


Table 3 



Here DIS means deep inelastic scattering, Bj stands for the Bjorken, and GLS for the Gross-Llewellyn Smith 
sum rules. The xF^ result is that of ref. 10. 

The previously existing average value, also taking into account NLO calculations, was 

as(M|) = 0.118 ±0.006; 

when including both our result and that of ref. 10 the new average and error become 

as(M|) = 0.1174 ± 0.0016. (5.2) 

We add a further comment. If we only kept, in the old determinations of a^, processes with spacelike 
momenta then, as noted by Bethke,[23] a slightly smaller value was obtained for the coupling: 

as{Mz) = 0.114 ± 0.005 (spacelike momenta). 

Our results correct this to some extent, by increasing slightly the average value for "spacelike" determinations. 

To finish this paper, we discuss briefly why our calculation yields such accurate results. First of all, 
and compared with neutrino data, our experimental input turns out to be much more precise, and also it 
extends over a much wider range of values. Actually, this wide range of values is one of the assets in our 
evaluation, particularly when comparing it to Z, r decays where only one value of is essentially available. 
Our range is such that not only leading logs, but next-to- leading logs vary appreciably: cf Table 1. The 
inclusion of NNL logs then stabilizes the results and decreases the errors and (very slightly) the x^/d.o.f.. 
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